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1.  Objectives 


There  were  four  main  objectives  in  the  project:  1)  develop  and  refine  classical  force  fields 
for  modeling  the  solid  and  liquid  phases  of  a  range  of  ionic  liquids;  2)  use  these  force 
fields  to  compute  thermodynamic  and  dynamic  properties  of  ionic  liquids  in  the 
condensed  phase;  3)  develop  new  simulation  methods  that  can  enable  property 
calculations  to  be  performed  (at  all  or  more  efficiently);  4)  systematically  study  how  the 
structure  and  composition  of  the  different  ionic  liquids  affects  the  thermodynamic  and 
transport  properties  of  ionic  liquids,  and  communicate  these  findings  to  synthetic  groups. 


2.  Status  of  Effort  -  Executive  Summary 

Objective  1:  New  force  fields  were  developed  for  the  following  molecular  compounds: 
hydrazine,  monomethylhydrazine  and  unsymmetrical  dimethyl  hydrazine.  Force  fields 
were  also  developed  for  the  following  ionic  liquids:  l-ethyl-3-methylimidaziolium 
ethylsulfate,  l-ethyl-3-methylimidazolium  butanesulfonate,  1 -ethyl-3 -methylimidazolium 
tricyanomethanide,  hydrizinium  nitrate,  methylhydrizinium  nitrate,  unsymmetrical 
dimethylhydrizinium  nitrate  and  2-hydroxyethylhydrizinium  nitrate.  As  part  of  a  study  on 
melting  points,  force  fields  were  developed  for  sodium  nitrate,  lithium  nitrate  and 
potassium  nitrate.  Finally,  a  set  of  scripting  tools  was  created  that  enables  the  generation 
of  new  ionic  liquid  force  fields  much  faster  than  was  possible  before. 

Objective  2:  Using  these  force  fields,  pure  liquid  properties  were  computed,  including 
density,  coefficient  of  thermal  expansion,  isothermal  compressibility,  cohesive  energy 
density,  vapor  pressure,  vapor-liquid  coexistence  curves,  critical  points,  and  fluid 
structure.  The  solubility  of  several  gases  in  ionic  liquids  were  computed,  as  were  the 
melting  points  and  crystal  polymorph  stability  of  one  particular  ionic  liquid. 

Objective  3:  A  new  Monte  Carlo  method  called  continuous  fractional  component  Monte 
Carlo  (CFC  MC)  was  developed  and  implemented  in  three  new  open  ensembles.  It  was 
used  to  compute  pure  and  mixed  gas  solubility  in  ionic  liquids.  The  pseudo-supercritical 
path  sampling  approach  was  refined  for  melting  point  and  polymorph  stability 
calculations.  It  was  applied  it  to  the  calculation  of  the  melting  point  of  an  ionic  liquid. 
Finally,  a  new  procedure  was  developed  for  computing  the  solubility  limit  of  crystalline 
materials  in  liquid  solvents.  It  was  tested  and  validated  on  sodium  chloride  solubility  in 
water  and  light  alcohols. 

Objective  4:  It  was  shown  how  hydrogen  bonding  between  water  and  the  anion  of 
different  ionic  liquids  slows  the  dynamics  and  can  contribute  to  higher  than  expected 
viscosities.  The  liquid  structure  of  the  energetic  ionic  liquid  2-hydroxyethylhydrizinium 
nitrate  was  examined.  It  was  shown  that  hydrogen  bonding  networks  form  that  result  in 
the  low  mobility  of  this  ionic  liquid.  A  combined  experimental  and  simulation  study  was 
carried  out  to  investigate  the  nature  of  ionic  liquids  in  the  vapor  phase.  It  was  determined 
that  the  most  stable  gas  phase  entities  of  ionic  liquids  are  ion  pairs.  Studies  to  examine 
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the  nature  of  the  products  that  form  when  ionic  liquids  are  subjected  to  strong  electric 
fields,  as  occurs  in  electrospray  experiments,  were  initiated. 

3.  Accomplishments  /  New  Findings 

3.1  Ionic  Liquids  in  the  Vapor  Phase 

A  series  of  simulations  were  carried  out  to  investigate  the  nature  of  ionic  liquids  in  the 
vapor  phase.  The  enthalpy  of  vaporization  of  different  sized  clusters  (both  neutral  and 
charged)  was  computed.  It  was  shown  that  by  far  the  most  energetically  favorable  species 
in  the  gas  phase  is  a  singe  ion  pair  that  forms  a  neutral  cluster  [1],  Simulations  in  support 
of  the  experimental  measurements  of  gas  phase  clusters  in  Prof.  Stephen  Leone’s  group 
were  also  carried  out[2].  Based  on  the  experiments  and  calculations,  it  was  conclusively 
shown  that  normal  vaporization  of  ionic  liquids  takes  place  by  neutral  clusters  of  a  single 
cation  and  anion  leaving  the  liquid  phase. 


Figure  1:  Calculated  and  experimental  enthalpy  of  vaporization  for  1 -alkyl-3  - 
methylimizadolium  bis(trifluoromethylsulfonyl)imide  ionic  liquids  as  a  function  of  the 
length  of  the  alkyl  group.  Despite  the  scatter  in  the  experimental  data,  the  simulations 
appear  to  do  a  good  job  capturing  both  the  magnitude  and  trend  in  the  data. 
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Figure  2:  Structure  of  the  ionic  liquid  l-ethyl-3-methylimidazolium 
bis(trifluoromethylsulfonyl)imide. 


Toward  the  end  of  the  project,  the  behavior  of  ionic  liquid 
droplets  in  vacuum  with  and  without  the  presence  of 
external  electric  fields  were  examined.  For  these 
simulations  the  ionic  liquid  l-ethyl-3-methylimidazolium 
bis(trifluoromethylsulfonyl)imide  ([emim][Tf2N])  was 
studied.  This  ionic  liquid  was  chosen  because  it  has  been 
studied  experimentally  a  great  deal  by  collaborators  Y. 
Chiu  and  co-workers  as  Hanscom  Air  Force  Base.  Figure  3 
shows  a  snapshot  from  one  of  these  simulations.  It  has  been 
found  that  for  this  ionic  liquid,  a  critical  field  strength  of 
about  1  V/nm  is  required  to  observe  emission,  in  good 
agreement  with  experimental  field  strengths. 


The  simulations  enabled  the  calculation  of  the  distribution 
of  ions  as  a  function  of  field  strength  and  droplet  size. 
Above  a  droplet  size  of  about  25  ion  pairs,  the  critical  field 
strength  for  emission  is  independent  of  droplet  size. 


Figure  3:  Droplet  of 
[emim] [Tf2N]  in  an  electric 
field.  Notice  the  emission  of  an 
ion  on  the  left  side  of  the 
droplet. 


3.2  Solubility  and  Dynamics  of  Ionic  Liquids  with  Water 

A  series  of  studies  were  conducted  to  compute  the  solubility  of  water  in  ionic  liquids  and 
to  see  how  water  affects  the  properties  of  ionic  liquids.  All  ionic  liquids  absorb  water  to 
some  extent,  and  so  it  is  critical  to  know  both  how  much  water  is  absorbed  and  how  the 
properties  are  altered  by  water.  Using  the  ionic  liquid  1 -ethyl-3 -methylimidazolium 
bis(trifhioromethylsulfonyl)imide  (shown  in  Figure  2)  the  isotherm  of  water  at  298  K  was 
computed.  Figure  4  shows  the  result.  The  only  experimental  data  available  was  the 
saturation  amount  of  water,  and  as  can  be  seen  the  simulation  results  are  consistent  with 
the  experiment.  To  perform  the  simulations,  a  new  Monte  Carlo  procedure  was  developed 
(see  section  3.3). 
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Figure  4:  Computed  water  isotherm  in  l-ethyl-3-methylimidazolium 

bis(trifluoromethylsulfonyl)imide.  The  black  open  symbol  is  the  only  experimental  data 
point  available.  From  ref.  [3J. 


It  has  been  claimed  previously  that  neglect  of  polarizability  leads  to  inaccuracies  in  the 
computed  transport  properties  of  ionic  liquids  such  as  viscosities  and  diffusivities.  No 
doubt  polarizability  does  yield  greater  mobility,  but  it  was  hypothesized  that  the 
consistently  low  viscosity  that  had  been  observed  by  many  previous  researchers  was  due 
more  to  the  way  in  which  they  computed  the  viscosity  than  to  anything  having  to  deal 
with  the  force  field.  That  is,  most  groups  rely  on  relatively  short  (100-500  ps)  simulations 
and  evaluate  the  viscosity  via  conventional  Green-Kubo  integration.  In  this  work  it  was 
shown  that  relevant  relaxation  times  are  much  longer  than  these  simulations,  and  so  the 
dynamics  being  computed  from  a  short  simulation  may  be  only  local  and  not  global. 
Thus,  a  proper  simulation  should  be  for  longer  than  the  longest  relevant  relaxation 
process.  To  overcome  this,  the  reverse  nonequilibrium  molecular  dynamics  (RNEMD) 
procedure  was  utilized  to  compute  the  viscosities  of  ionic  liquids.  Figure  5  shows  how 
are  calculations  match  experiments.  Note  that,  even  with  a  fixed  charge  model,  the 
agreement  between  experiment  and  simulation  is  excellent. 
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Figure  5:  Computed  viscosity  as  a  function  of  temperature  for  l-ethyl-3- 
methylimidazolium  bis(trifluoromethylsulfonyl)imide.  The  black  line  is  a  VFT  fit  to  a 
large  number  of  experimental  data  [3], 

Self-diffusivities  were  also  computed,  and  these  are  shown  in  Figure  6.  The  computed 
dynamics  are  about  a  factor  of  two  slower  than  the  NMR  results  of  Vasenkov  et  ah, 
which  is  consistent  with  the  slightly  higher  computed  viscosity  observed  in  Figure  5.  The 
trend  with  temperature  is  captured  well,  as  is  the  curious  result  that  the  larger  cation  has  a 
larger  self-diffusivity  than  the  smaller  anion. 
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Figure  6:  Computed  self-diffusivity  for  l-ethyl-3-methylimidazolium  cation  and 
bis(trifluoromethylsulfonyl)imide  anion  as  a  function  of  temperature.  From  ref  [13]. 

3.3  New  Simulation  Methods 

Three  new  simulation  methods  were  developed  and  applied  to  compute  properties  of 
ionic  liquids.  The  first  method,  called  continuous  fractional  component  Monte  Carlo 
(CFC  MC)  enables  efficient  transfer  of  molecules  between  systems  in  open  ensemble 
Monte  Carlo.  The  method  was  described  in  references  4  and  8,  and  applied  to  compute 
CO2  solubility  in  an  ionic  liquid  in  reference  7.  Figure  4  shows  water  solubility  results 
computed  with  this  method,  while  Figure  7  shows  a  comparison  between  calculated  and 
experimental  isotherms  for  CO2  and  l-n-hexyl-3-methylimidazolium 
bis(trifluoromethylsulfonyl)imide.  Agreement  is  quantitative  in  both  cases. 
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Figure  7:  Calculated  (red  circles)  and  experimental  (open  symbols)  solubility  of  C02  in 
the  ionic  liquid  1  -n-hexyl-3-methylimidazolium  bis(trifluoromethylsulfonyl)imide  at  323 
K.  From  ref.  [7] . 


The  method  has  been  implemented  in  the  osmotic  and  grand  ensembles  [4]  as  well  as  the 
Gibbs  ensemble  [8].  The  method  is  being  implemented  in  an  open-source  Monte  Carlo 
code  in  order  to  make  it  more  widely  available  to  the  broader  research  community. 

The  second  new  simulation  method  of  this  project  was  a  refinement  of  an  earlier 
technique  called  pseudosupercritcal  path  sampling  (PSPS).  This  method  enables  one  to 
compute  melting  points  and  polymorph  stability.  The  method  works  as  follows.  Relative 
free  energy  curves  for  a  crystal  and  the  liquid  phase  of  the  same  substance  are  computed 
by  conducting  isotherm-isobaric  simulations  of  both  phases.  The  temperature-dependence 
of  the  enthalpy  is  computed  and  the  relative  free  energies  are  obtained  using  the  Gibbs- 
Helmholtz  equation.  Note  that  this  gives  the  free  energy  of  the  solid  and  liquid  phases 
relative  to  some  unknown  integration  constant.  As  an  example,  Figure  8  shows  such  a 
construction  for  the  ionic  liquid  l-n-butyl-3-methylimidazolium  chloride. 

The  next  step  is  to  determine  the  constants  for  the  two  phases  and  in  so  doing,  obtain  the 
free  energies  of  both  phases  as  a  function  of  temperature.  When  the  two  are  equal,  the 
system  is  in  coexistence.  To  do  this,  a  thermodynamic  cycle  is  constructed  between  the 
two  phases  at  a  given  temperature. 
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Figure  8:  (left)  -H/RT2  for  the  liquid  and  crystal  phases  of  l-n-butyl-3- 
methylimidazolium  chloride,  (right)  The  solid  and  liquid  free  energies  using  a  reference 
state  at  300K. 


A  schematic  of  the  cycle  is  shown  in  Figure  9. 


State  2:  Weak  fluid 


State  4:  Ordered 
weak  phase 


Figure  9:  Thermodynamic  cycle  for  determining  the  free  energy  difference  between  a 
crystal  and  a  liquid  at  a  given  temperature. 


In  going  from  state  1  to  state  2,  a  liquid  is  taken  to  a  weakly  interacting  liquid  by  scaling 
the  intermolecular  energy  between  atoms.  From  state  2  to  state  3,  the  density  of  the 
weakly  interacting  fluid  is  adjusted  to  yield  the  crystal  density.  State  4  is  formed  by  using 
tethering  potentials  that  order  the  fluid  into  the  position  found  in  a  crystal.  Finally,  state  5 
(the  crystal  phase)  is  recovered  by  turning  back  on  the  intermolecular  interactions  and 
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removing  the  tethering  potential.  At  each  step  the  free  energy  change  is  computed  via 
thermodynamic  integration.  The  result  is  the  actual  free  energy  difference  between  the 
liquid  and  crystal  at  a  particular  temperature.  Given  this  difference  at  one  temperature, 
the  difference  between  the  integration  constants  for  the  two  phases  are  known,  and  the 
curves  in  Figure  8  can  be  shifted  to  find  the  coexistence  temperature. 

The  PSPS  method  was  used  to  conduct  the  first  rigorous,  free  energy-based  calculation  of 
the  melting  point  of  an  ionic  liquid,  and  the  first  calculation  in  which  the  thermodynamic 
stability  of  two  ionic  liquid  crystal  polymorphs  were  computed  [6].  Figure  10  shows  two 
different  polymorphs  of  l-n-butyl-3-methylidazolium  chloride. 


Figure  10:  Orthorhombic  (left)  and  monoclinic  (middle)  structures  of  the  l-n-butyl-3- 
methylidazolium  chloride  ([bmimJCl)  ionic  liquid.  The  bmim  cation  is  shown  on  the 
right.  The  only  difference  between  these  two  ionic  liquids  is  the  torsion  angle  about  the 
C7  and  C8  carbons. 

The  table  below  shows  the  results  of  the  calculations. 


polymorph 

T~pt(K) 

T£Lr(K) 

AHr*(kJ/mol) 

AHfff  (kJ/mol) 

AHf“JJ(kJ/mol) 

orthorhombic 

337-339 

302  ±5 

308  ±9 

18.5-21.5 

26  ±3 

19  =4 

monoclinic 

318-340 

NA 

333  ±6 

9.3-14.5 

NA 

23  =3 

The  “CLP”  results  are  those  obtained  using  a  force  field  from  Canongia-Lopes  and 
Padua,  while  the  “CMJJ”  is  the  potential  developed  in  this  work.  The  CLP  force  field  did 
not  yield  a  stable  monoclinic  crystal.  The  CMJJ  force  field  does  a  good  job  matching  the 
melting  points  and  heats  of  fusion. 
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The  method  was  then  applied  to  compute  the  melting  points  of  LiNOs,  NaN03  and  KNO3 
[17].  These  systems  were  chosen  because  there  experimental  data  for  them  and  they 
make  a  good  system  to  test  the  method  on.  Figure  1 1  shows  the  free  energy  difference 
between  the  crystal  phase  and  liquid  phase  as  a  function  of  temperature  for  LiNC>3.  The 
computed  melting  point  is  459  K  while  the  experimental  value  is  531  K.  The  difference  is 
due  to  inaccuracies  in  the  force  field,  though  it  is  clear  from  how  flat  the  curve  is  in 
Figure  11  that  differences  in  free  energy  of  a  few  kJ/mol  can  result  in  differences  in 
melting  points  of  10s  of  degrees.  This  points  out  how  difficult  it  is  to  get  very  accurate 
melting  points. 


Figure  11:  Computed  free  energy  difference  between  the  crystal  and  liquid  phases  of 
LiNOs. 

The  third  method  developed  in  this  project  builds  on  the  PSPS  method  and  was  used  to 
compute  the  solubility  limit  of  crystals  in  liquids  [18].  Figure  12  shows  a  schematic  of 
how  the  method  works,  where  species  1  is  the  crystalline  solute  and  species  2  is  the  liquid 
solvent.  In  this  case,  the  PSPS  method  is  used  to  compute  the  difference  in  chemical 
potential  between  the  crystalline  phase  and  a  weakly  interacting  fluid  phase  (Ajfh). 
Expanded  ensemble  methods  are  used  to  compute  the  chemical  potential  difference 
between  the  crystalline  species  1  in  a  common  reference  state  (typically  an  ideal  gas)  and 
the  weak  fluid  (A^i).  In  addition,  the  chemical  potential  difference  (A[f))  is  computed, 
which  corresponds  to  the  chemical  potential  difference  between  the  solute  in  the  common 
reference  state  and  in  solution  at  a  given  concentration.  Figure  13  shows  a  typical  result 
for  sodium  chloride  in  water.  Note  the  huge  free  energy  difference  between  the  reference 
state  and  the  solution  phase.  This  gives  an  indication  of  how  difficult  it  is  to  determine 
free  energies  of  salts  in  aqueous  solutions.  The  only  way  this  calculation  was  possible 
was  to  break  up  the  simulation  into  five  “windows”  where  different  parts  of  the  potential 
were  “turned  on”  at  different  times.  In  addition,  advanced  biasing  methods  were  used  to 
improve  the  transitions  between  adjacent  states  in  the  expanded  ensemble  [18]. 
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By  computing  this  quantity  for  a  series  of  different  solute  concentrations,  one  can  close 
the  thermodynamic  cycle  in  Figure  12  and  determine  the  concentration  that  corresponds 
to  coexistence.  The  temperature  dependence  can  be  incorporated  by  repeating  the 
calculation  at  different  temperatures.  Figure  13  shows  an  example  of  one  such 
calculation. 


Figure  12:  A  schematic  of  the  method  to  determine  the  solubility  limit  of  crystalline 
materials.  The  chemical  potential  of  both  phases  is  determined  relative  to  a  common 
reference  state,  such  that  relative  (and  not  absolute)  chemical  potentials  are  required. 


Figure  13:  Excess  chemical  potential  of  sodium  chloride  in  water  at  0.2 1M  as  a  function 
of  the  coupling  parameter.  Different  line  colors  and  types  are  used  to  distinguish 
between  the  results  of  the  five  sampling  windows  used. 


3.4  Simulations  of  Energetic  Ionic  Liquids 

A  wide  range  of  simulations  on  ionic  liquids  based  on  the  hydrizinium  cation  paired  with 
the  nitrate  anion  were  carried  out  [11].  Hydrizinium  nitrate  ([N2H5][N03]), 
monomethylhydrizinium  nitrate  ([CH3N2H4][N03]),  and  unsymmetrical 
dimethylhydrizinium  nitrate  ([(CH3)2N2H3][N03])  have  melting  points  below  100°C,  and 
thus  may  be  considered  as  ionic  liquids. 

To  start  with,  force  fields  for  molecular  hydrazine-like  compounds  were  developed.  The 
phase  behavior  of  these  compounds  was  computed  using  Gibbs  ensemble  Monte  Carlo. 


11 


Figure  14  shows  the  computed  and  experimental  vapor-liquid  coexistence  curve  for 
hydrazine  (left)  and  enthalpy  of  vaporization  (right).  Agreement  with  experiment  is 
excellent. 


700 

*  600 

o 

3  500  • 

♦ 

^  400  ♦ 
O  300  # 


♦  sim  N=31Q  exp  sim  N=51 0 1 


200 


0  0  0  2  0  4  0  6  0  8  1  0  1  2 

density  [g/cm3] 


50 

^  45 
O  40 

^  35 

“  30 

20 

15 


-I  lilrh 

♦  sim 

f'f  | 

P 

cj  exp 

1 

$ 

200  300  400  500  600 

T  [K] 


Figure  14:  Computed  and  experimental  vapor-liquid  coexistence  curve  for  hydrazine 
(left)  and  the  associated  enthalpy  of  vaporization  as  a  function  of  temperature  (right). 
Simulation  results  are  in  quantitative  agreement  with  experiments. 

The  Lennard-Jones  force  field  parameters  for  these  molecular  species  were  then 
transferred  to  the  ionic  liquid  analogues,  and  a  range  of  properties  were  computed  for  the 
ionic  liquids.  Initially,  simulations  of  the  crystal  phases  of  the  ionic  liquids  were  carried 
out.  Figure  15  shows  the  results. 
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Figure  15:  Simulated  crystal  phases  of  (from  left  to  right):  hydrizinium  notrate, 
monomethylhydrizinium  nitrate  and  unsymmetrical  dimethylhydrizinium  nitrate. 


The  simulated  monoclinic  [N2H5][N03]  (Figure  15,  left)  shows  moderate  distortion 
compared  to  the  x-ray  lattice  structure.  Although  the  zig-zag  pattern  along  the  a-axis 
direction  is  maintained,  there  are  considerable  thermal  fluctuations.  The  lattice  contracts 
in  the  a  and  c  directions  by  5.5  and  2.9%,  respectively,  and  contracts  in  the  b  direction  by 
3.5%.  The  (3  angle  is  in  very  good  agreement  with  the  experimental  angle,  differing  by 
only  0.04%.  Overall,  despite  the  distortion  in  the  cell,  the  simulated  density  is  only  5.2% 
greater  than  the  experimental  value  and  agrees  reasonably  well. 


The  simulated  lattice  for  [CTfNifEjfNCf]  (Figure  15,  middle)  shows  major  distortions 
compared  to  the  experimental  x-ray  structure.  Thus,  the  force  field  is  unsuccessful  in 
simulating  the  stable  monoclinic  experimental  crystal  phase.  The  regular  repeating 
pattern  of  the  experimental  structure  is  no  longer  visible,  and  thermal  fluctuations  are 
dramatic.  As  seen  in  the  figure,  the  simulated  crystal  expands  in  the  a  direction  by  15.8% 
and  contracts  in  the  b  direction  by  11.1%.  Distortion  in  the  perpendicular  c  direction  is 
minor  at  2.0%.  Even  the  (3  angle  increases  by  over  1°  (1.3%).  The  nearly  equal  and 
opposite  distortions  in  a  and  b  axes  allow  the  simulated  density  to  remain  in  good 
agreement  with  the  experimental  density,  being  only  4.5%  lower,  an  entirely  fortuitous 
result.  In  the  experimental  structure  from  Figure  15,  considerable  free  volume  is 
observed,  and  upon  closer  inspection,  the  nitrate  anions  are  stacked  in  eclipsed  columns. 
However,  upon  simulation,  the  nitrates  remain  stacked,  but  in  an  irregular  fashion, 
accompanied  by  a  decrease  in  the  free  volume  of  the  cell. 
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The  simulated  lattice  for  [(CFf^^FbJfNOs]  (Figure  15,  right)  is  in  excellent  agreement 
with  the  experimental  orthorhombic  x-ray  structure.  Aside  from  small  thermal 
fluctuations,  the  simulated  lattice  retains  its  shape  and  symmetry  relative  to  the  crystal 
lattice.  Deviations  in  the  a,  b,  and  c  axes  are  are  +1.3%,  +1.4%,  and  -3.6%,  respectively. 
The  cell  angles  were  not  allowed  to  vary  from  their  90°  ideal  values.  The  simulated 
density  is  only  1.0%  lower  than  the  experimental  density. 


Interestingly,  the  nitrate  anions  in  the  experimental  lattices  of  [N2H5][N03]  and 
[(CH3)2N2H3][N03]  are  not  stacked  in  an  apparent  regular  fashion  as  in  [CFf^FEJfNOs], 
but  rather  appear  staggered.  The  simulated  lattices  of  the  former  two  crystals, 
coincidentally  or  not,  show  the  smallest  deviations  from  experiment,  although  in 
[N2H5][N03]  some  may  argue  that  the  deviations  are  important.  Nonetheless,  a  possible 
reason  for  the  unstable  [CH3N2H4][N03]  simulated  phase  may  be  due  to  limitations  in  the 
nitrate  force  field,  and  the  driving  force  to  maximize  hydrogen  bonding  interactions  via 
irregularities  in  the  stacking  of  the  anion.  Otherwise,  the  inability  of  the 
methylhydrazinium  force  field  to  perform  well,  considering  the  exceptional  performance 
for  the  analogous  neutral  MEH,  is  inexplicable. 

Figure  16  shows  computed  liquid  densities  and  molar  volumes  for  the  three  hydrizinium 
salts. 
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Figure  16:  Computed  liquid  densities  and  molar  volumes  for  the  hydrizinium  salts,  in 
addition  to  snapshots  showing  typical  interactions  between  the  cations  and  anions. 


The  liquids  were  simulated  at  348  K,  above  the  melting  points  of  all  three  salts.  All  three 
have  very  high  densities  at  the  simulated  temperature,  decreasing  from  1.681  g/cm3  for 
[N2H5][N03]  to  1.230  g/cm3  for  [(CH3)2N2H3][N03].  As  the  density  decreases,  the  molar 
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volume  increases  down  the  hydrazinium  series.  The  molar  volume  of 
[(CH3)2N2H3][N03]  is  nearly  twice  that  of  [N2H5][N03].  As  shown  below,  the  dynamics 
of  these  salts  is  sluggish,  and  we  believe  this  is  due  in  part  to  the  strong  hydrogen 
bonding  that  takes  place  between  cations  and  the  nitrate  anions.  Figure  16  shows 
snapshots  from  the  simulations  where  such  hydrogen  bonding  interactions  can  be 
observed. 

The  self-diffusion  coefficients  (Ds)  and  rotational  time  constants  were  also  evaluated  for 
the  cation  and  anion  of  each  salt.  In  each  salt,  the  anion  exhibits  faster  translations 
dynamics  that  the  cation  (larger  Ds  values).  Additionally,  [N2H5][N03]  exhibits  faster 
dynamics  overall  relative  to  [CH3N2H4][N03],  and  [(CH3)2N2H3][N03],  and  all  are  on  the 
order  of  10'  m  /s,  although  the  3  ns  simulations  may  be  inadequate  to  sample  the 
diffusive  regime,  given  the  slow  dynamics.  Rotational  time  constants  for  the  cations  and 
anions  of  each  IL  are  on  the  order  of  several  hundred  picoseconds. 

The  hydrazinium  salt,  2-hydroxyethylhydrazinium  nitrate  (2-HEHN)  is  a  low  melting 
ionic  liquid  that  has  considerable  potential  for  energetic  applications.  Although  little 
experimental  data  is  available  in  the  literature,  predictions  were  made  concerning  the 
physical  and  dynamical  properties  of  this  IL  based  on  previous  validation  of  the 
molecular  force  field.  As  shown  in  Figure  17,  the  liquid  densities  decrease  linearly  with 
increasing  temperature.  The  experimental  density  at  298  K  of  1.42  g/cm  (Brand,  A.  J.; 
Drake,  G.W.  “Energetic  hydrazinium  salts.”  U.S.  patent  #6218577)  is  in  excellent 
agreement  with  the  simulated  value  of  1.428  g/cm  .  The  liquid  structure  was  evaluated 
via  radial  distribution  analysis  of  the  system  components.  The  cation-cation,  cation- 
anion,  and  anion- anion  COM-COM  radial  distribution  functions  are  shown  in  Figure  18. 
The  attractive  cation-anion  (+-)  g(r)  is  the  most  intense  and  occurs  at  the  shortest  distance 
of  4  A.  A  significant  cation-anion  longer  ranger  interaction  occurs  at  ~9  A.  The 
repulsive  cation-cation  (++)  and  anion-anion  (— )  interactions  are  less  intense  and  occur  at 
longer  distances  of  5.7  A  and  7.3  A,  respectively.  The  two  anion-anion  peaks  at  shorter 
distance  heave  intensity  less  than  1  and  are  not  important.  The  cation-cation  and  anion- 
anion  also  show  weaker  correlations  out  at  about  12  A. 
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Figure  17:  Computed  density  of  hydroxyethylhydrizinium  nitrate  as  a  function  of 
temperature. 
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Figure  18:  COM-COM  radial  distribution  functions  for  cation/cation  (++),  cation/anion 
(+-),  and  anion/anion  (—)  interactions  in  2-HEHN. 
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The  IL  2-HEHN  is  expected  to  exhibit  strong  hydrogen  bonding  in  the  liquid  phase.  As 
shown  in  Figure  19,  the  major  hydrogen-bonding  interactions  between  the  cation  and 
anion  have  been  evaluated  using  a  g(r)  correlation.  Here,  the  interaction  between  the 
nitrate  oxygen  atoms  and  either  the  hydroxyl  hydrogen  or  amine  hydrogens  of  2-HEH+ 
has  been  studied.  The  N-H  -O  interaction  is  far  weaker  than  the  O-H-O  interaction. 
The  former  occurs  at  ~1.9  A  while  the  latter  occurs  at  the  much  shorter  distance  of  ~1.5 
A.  The  1.5  A  hydrogen  bond  is  very  short  relative  to  hydrogen  bonds  in  general,  and 
suggests  that  this  interaction  is  very  strong  in  the  liquid.  This  is  supported  by  the 
sharpness  and  intensity  of  the  g(r)  curve  for  this  peak  relative  to  the  N-H  -  •  0  interaction. 
The  correlations  at  about  3.5  to  4  A  are  associated  with  the  neighboring  oxygens  of  the 
nitrate  anion.  The  information  provided  by  these  g(r)  curves  is  evidence  of  an  IL  that  is 
dominated  by  hydrogen  bonding  due  to  the  presence  of  many  significant  donor/acceptor 
sites. 


r(A) 


Figure  19:  Site-site  radial  distribution  functions  for  hydroxyl:  nitrate  interactions  and 
amine:nitrate  interactions  in  2-HEHN. 


The  dynamics  of  2-HEHN  were  probed  by  calculating  the  MSD  vs.  time  for  the  cation 
and  anion,  as  well  as  their  self-diffusion  coefficients.  The  temperature  dependence  of  the 
dynamics  was  evaluated  over  the  temperature  range  of  223  to  473  K.  Figure  20  shows 
the  MSD  (A2)  vs.  temperature  (K)  at  four  different  observation  times.  Information  about 
the  melting/glass-transition  region  for  an  IL  can  be  elucidated  from  the  MSD  vs. 
temperature  curve  in  the  region  where  the  ions  exhibit  small  fluctuations  (<0.1  A2), 
indicating  lack  of  global  motion.  In  Figure  20,  the  MSDs  increase  rapidly  at  higher 
temperature,  making  it  difficult  to  establish  behavioral  characteristics  at  lower 
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MSD  (A2) 


temperatures.  Thus,  the  breakout  plots  illustrate  behavior  over  the  entire  temperature 
range  focusing  on  smaller  MSDs.  As  seen  in  the  two  figures  for  the  cation  and  anion 
below  about  250  K,  motion  over  the  1000  ps  observation  range  is  <0.1  A2,  indicating 
glassy  behavior.  However,  above  250  K,  displacements  over  the  1000  ps  timeframe  are 
greater  than  0. 1  A  and  mark  the  start  of  departure  from  linear  behavior.  Qualitatively, 
this  is  good  agreement  with  the  reported  melting  point  of  2-HEHN  of  <  248  K.  Further 
studies  of  melting  point  using  the  PSPS  method  will  allow  quantitative  confirmation. 


200  250  300  350  400  450 

T  (K) 


200  250  300  350  400  450 

T  (K) 


Figure  20:  Simulated  temperature  dependence  of  2-HEHN  cation  and  anion  MSDs  (A2)  at  several 


elavsed  times. 


Self-diffusion  coefficients  were  estimated  from  the  MSD  vs.  time;  at  the  lower 
temperatures,  the  cation  and  anion  of  2-HEHN  exhibit  sub-diffusive  behavior,  making  it 
impossible  to  quantitatively  assign  Ds.  At  low  temperatures  in  the  range  of  223  to  348  K, 
the  translational  motion  of  the  ions  is  very  small,  on  the  order  of  10'  m  Is,  much  slower 
than  previously  reported  ILs.  This  slow  behavior  is  likely  attributable  to  the  strong 
hydrogen-bonding  behavior  of  the  IL.  Although  it  is  difficult  to  appraise  these  results  in 
light  of  the  lack  of  experimental  evidence,  it  seems  plausible  that  the  slow  dynamics  may 
be  due  to  hydrogen  bonds  that  are  too  strong.  At  temperatures  of  373  K  and  above,  a 
speed-up  in  the  dynamics  is  observed.  At  the  highest  temperature  of  473  K,  the  MSDs 
are  ~4  xlO'11  m2/s  and  ~6xl  ()"' 1  m2/s  for  the  cation  and  anion,  respectively.  The  hydrogen 
bonds  are  labile  this  temperature  region,  allowing  the  species  to  move  more  freely  in  the 
bulk.  Employing  dampened  charges  or  using  a  polarizable  may  aid  in  the  future  study  of 
these  materials. 
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